load("simulation_output/main_sims.Rdata", .GlobalEnv)

summary_sims <- as.data.frame(do.call(rbind, lapply(1:length(sims_m_params), FUN = function(i) sims_m_params[[i]]$summary_sims)))

color_vec <- c('orange', 'orange3', 'dodgerblue3', 'navyblue', 'indianred3')

p1 <- ggplot(summary_sims) + 
   geom_hline(aes(yintercept = 3), linetype = 'dashed', col = 'grey59') + 
  geom_point(aes(x = S1, y = b1_est - 12), col = color_vec[1]) +
geom_line(aes(x = S1, y = (b1_est - 12)), col = color_vec[1]) +
geom_point(aes(x = S1, y = (b1_est_tilde - 12)), col = color_vec[3]) +
geom_line(aes(x = S1, y = (b1_est_tilde - 12)), col = color_vec[3]) +
geom_point(aes(x = S1, y = b2_est + 3), col = color_vec[2]) +
geom_line(aes(x = S1, y = (b2_est + 3)), col = color_vec[2]) +
geom_point(aes(x = S1, y = (b2_est_tilde + 3)), col = color_vec[4]) +
geom_line(aes(x = S1, y = (b2_est_tilde + 3)), col = color_vec[4])  +
labs(y = 'bias', x = TeX('$\\S_1$')) + 
  annotate("text", x = 3.8, y = -9, label = TeX('$\\b_1$'), colour = color_vec[1]) + 
  annotate("text", x = 0.5, y = 2, label = TeX('$\\b_2$'), colour = color_vec[2]) + 
  theme_bw() +
  theme(panel.grid.major = element_blank()) +
  annotate("text", x = 3.85, y = 0.8, label = TeX('$\\tilde{\\beta}_1$'), colour = color_vec[3]) + 
  annotate("text", x = 3.8, y = -0.8, label = TeX('$\\tilde{\\beta}_2$'), colour = color_vec[4]) + 
  labs(y = 'Bias', y = TeX('$\\S_{1}$')) + ylim(c(-12, 5)) + 
  annotate("text", x = 1, y = 3.9, label = paste0('Wrong sign on b2'), col = 'grey59')


p2 <- ggplot(summary_sims) + geom_point(aes(x = S1, y = b1_se_est_tilde), col = color_vec[4]) +
  geom_line(aes(x = S1, y = b1_se_est_tilde), col = color_vec[4]) +
  geom_line(aes(x = S1, y = b1_se_true_tilde), col = color_vec[3]) +
    geom_point(aes(x = S1, y = b1_se_true_tilde), col = color_vec[3]) + 
  geom_point(aes(x = S1, y = b2_se_est_tilde), col = color_vec[4]) +
  geom_line(aes(x = S1, y = b2_se_est_tilde), col = color_vec[4]) +
  geom_line(aes(x = S1, y = b2_se_true_tilde), col = color_vec[3]) +
  geom_point(aes(x = S1, y = b2_se_true_tilde), col = color_vec[3]) +
  theme_bw() + 
  theme(panel.grid.major = element_blank()) +
    annotate("text", x = 3.52, y = 0.48, label = TeX('$\\hat{SE}(\\beta_1)$'), colour = color_vec[4]) + 
    annotate("text", x = 3.67, y = 0.37, label = TeX('$\\SE(\\beta_1)$'), colour = color_vec[3]) +
 annotate("text", x = 3.6, y = 0.22, label = TeX('$\\hat{SE}(\\beta_2)$'), colour = color_vec[4]) + 
  annotate("text", x = 3.65, y = 0.12, label = TeX('$\\SE(\\beta_2)$'), colour = color_vec[3]) + 
  labs(y = 'Standard Error', x = TeX('$\\S_1$'))



summary_sims$MSE_b1 <- sqrt(summary_sims$b1_se_true^2 + (12 - summary_sims$b1_est)^2)
summary_sims$MSE_beta1 <- sqrt(summary_sims$b1_se_true_tilde^2 + (12 - summary_sims$b1_est_tilde)^2)
summary_sims$MSE_b2 <- sqrt(summary_sims$b2_se_true^2 + (-3 - summary_sims$b2_est)^2)
summary_sims$MSE_beta2 <- sqrt(summary_sims$b2_se_true_tilde^2 + (-3 - summary_sims$b2_est_tilde)^2)



p3 <-  ggplot(summary_sims) + 
geom_point(aes(x = S1, y = MSE_beta1), col = color_vec[3]) +
  geom_line(aes(x = S1, y = MSE_beta1), col = color_vec[3])  +
  geom_line(aes(x = S1, y = MSE_b1), col = color_vec[1]) +
  geom_point(aes(x = S1, y = MSE_b1), col = color_vec[1]) +
  geom_point(aes(x = S1, y = MSE_beta2), col = color_vec[4]) +
  geom_line(aes(x = S1, y = MSE_beta2), col = color_vec[4])  +
  geom_line(aes(x = S1, y = MSE_b2), col = color_vec[2]) +
  geom_point(aes(x = S1, y = MSE_b2), col = color_vec[2])  +
  annotate("text", x =  3.2, y = 9.8, label = TeX('$\\b_1$'), colour = color_vec[1]) + 
  annotate("text", x = 3.2, y = 1, label = TeX('$\\tilde{\\beta}_1$'), colour = color_vec[3]) +
  annotate("text", x =  3.8, y = 4.5, label = TeX('$\\b_2$'), colour = color_vec[2]) + 
  annotate("text", x = 3.8, y = -.5, label = TeX('$\\tilde{\\beta}_2$'), colour = color_vec[4]) +
  theme_bw() + 
  labs(x = TeX('$\\S_1$')) +
  theme(panel.grid.major = element_blank()) +
  labs(y = "Root Mean Squared Error") + geom_hline(aes(yintercept = 0), linetype = 'dashed', col = 'grey')



 load("simulation_output/time_test_analytical.Rdata", .GlobalEnv)

 p4 <- ggplot(time_variance) + 
   geom_point(aes(x = N, y = simulation, col = 'Simulation')) +
   geom_line(aes(x = N, y = simulation, col = 'Simulation')) +
   geom_point(aes(x = N, y = analytical, col = 'Analytical')) +
   geom_line(aes(x = N, y = analytical, col = 'Analytical')) +
   scale_color_manual(values = c('navyblue', 'orange')) + theme_bw() +
   labs(x = 'n', y = 'Compute Time (seconds)') +
   annotate("text", x = 4500000, y = 10, label = 'Simulation', colour = color_vec[2]) +
   annotate("text", x = 4000000, y = 60, label = 'Analytical', colour = color_vec[4]) +
   theme(legend.position="none") +
   theme(panel.grid.major = element_blank()) +
   scale_x_continuous(breaks = c(100000, 1000000, 2000000, 3500000, 5000000)) +
   geom_hline(aes(yintercept = 0), linetype = 'dashed', col = 'grey')



upperpanel <- suppressWarnings(ggarrange(p1, p3))

ggsave(upperpanel, file = 'plots/figure1.pdf', width = 9, height = 4)

lowerpanel <- suppressWarnings(ggarrange(p2, p4))

ggsave(lowerpanel, file = 'plots/figure2.pdf', width = 9, height = 4)


